In [1]:
%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt

In [2]:
# implementation of pairwise distance function

# make datasets of coordinates
I = -np.random.randint(10,size=(4,2))
J = -np.random.randint(10,size=(3,2))

I


Out[2]:
array([[-4, -5],
       [-9, -3],
       [-7, -5],
       [-2, -5]])

In [3]:
J


Out[3]:
array([[-4, -2],
       [ 0, -2],
       [-6, -8]])

In [4]:
from scipy.spatial import distance

distance.cdist(I, J)


Out[4]:
array([[ 3.        ,  5.        ,  3.60555128],
       [ 5.09901951,  9.05538514,  5.83095189],
       [ 4.24264069,  7.61577311,  3.16227766],
       [ 3.60555128,  3.60555128,  5.        ]])

In [23]:
# ignoring other metrics, we want to implement euclidean distance
# sqrt((x0-x1)**2 + (y0-y1)**2)

In [12]:
np.mgrid[:I.shape[0],:J.shape[0]]


Out[12]:
array([[[0, 0, 0],
        [1, 1, 1],
        [2, 2, 2],
        [3, 3, 3]],

       [[0, 1, 2],
        [0, 1, 2],
        [0, 1, 2],
        [0, 1, 2]]])

In [20]:
# create a rectangular grid that is Ni by Nj
i, j = np.mgrid[:I.shape[0],:J.shape[0]]
i


Out[20]:
array([[0, 0, 0],
       [1, 1, 1],
       [2, 2, 2],
       [3, 3, 3]])

In [25]:
j


Out[25]:
array([[0, 1, 2],
       [0, 1, 2],
       [0, 1, 2],
       [0, 1, 2]])

In [24]:
# now we want the pairs of points I[i] and J[j]
J[j]


Out[24]:
array([[[-4, -2],
        [ 0, -2],
        [-6, -8]],

       [[-4, -2],
        [ 0, -2],
        [-6, -8]],

       [[-4, -2],
        [ 0, -2],
        [-6, -8]],

       [[-4, -2],
        [ 0, -2],
        [-6, -8]]])

In [32]:
J[j]


Out[32]:
array([[[-6, -8],
        [-1, -7],
        [-7, -2]],

       [[-6, -8],
        [-1, -7],
        [-7, -2]],

       [[-6, -8],
        [-1, -7],
        [-7, -2]],

       [[-6, -8],
        [-1, -7],
        [-7, -2]]])

In [25]:
# to subtract x0-x1 and y0-y1 simultaneously, we need this expression
I[i] - J[j]


Out[25]:
array([[[ 0, -3],
        [-4, -3],
        [ 2,  3]],

       [[-5, -1],
        [-9, -1],
        [-3,  5]],

       [[-3, -3],
        [-7, -3],
        [-1,  3]],

       [[ 2, -3],
        [-2, -3],
        [ 4,  3]]])

In [28]:
# then we square it
s = (I[i] - J[j])**2
s


Out[28]:
array([[[ 0,  9],
        [16,  9],
        [ 4,  9]],

       [[25,  1],
        [81,  1],
        [ 9, 25]],

       [[ 9,  9],
        [49,  9],
        [ 1,  9]],

       [[ 4,  9],
        [ 4,  9],
        [16,  9]]])

In [31]:
# now take the sum of squares along the 3rd axis (we're summing (x0-x1)**2 and (y0-y1)**2)
sq = np.sum(s,axis=2)
sq


Out[31]:
array([[ 9, 25, 13],
       [26, 82, 34],
       [18, 58, 10],
       [13, 13, 25]])

In [32]:
# now take the square root
np.sqrt(sq)


Out[32]:
array([[ 3.        ,  5.        ,  3.60555128],
       [ 5.09901951,  9.05538514,  5.83095189],
       [ 4.24264069,  7.61577311,  3.16227766],
       [ 3.60555128,  3.60555128,  5.        ]])

In [33]:
# here it is as a function
def cdist(I, J):
    i, j = np.mgrid[:I.shape[0],:J.shape[0]]
    return np.sqrt(np.sum((I[i] - J[j])**2,axis=2))

cdist(I, J)


Out[33]:
array([[ 3.        ,  5.        ,  3.60555128],
       [ 5.09901951,  9.05538514,  5.83095189],
       [ 4.24264069,  7.61577311,  3.16227766],
       [ 3.60555128,  3.60555128,  5.        ]])